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The  multiconf igurational  density  functional  theory  (MCDFT)  software  Schrodinger,  Inc. 
has  developed  in  Phase  I  of  this  proejct  combines  many  of  the  advantages  of  gener¬ 
alized  valence  bond  (GVB)  and  restricted  configuration  interaction  (RCI)  techniques 
with  those  of  DFT.  WE  have  explored  several  avenues  for  improving  the  implementa¬ 
tion  of  MCDFT  mehtods  and  ahve  incorporated  the  MCDFT  code  into  a  development  version 
of  our  commercially  successful!  electronic  structure  program,  PS-CVB  [1],  The  work  we 
have  performed  in  Phase  I  has  centered  upon  these  five  tasks;  partitioning  of  the 
CVB  two-electron  energy  into  Coulomb,  exchange,  and  intra-pair  terms;  coding  and 
developing  CVB-RCI-DFT  (post-SCF  DFT  on  a  CVB-RCI  density);  inclusion  of  self-inter- 
action-corrected  (SIC)  density  functionals;  optimization  of  hybrid  MCDFT  methods; 
and  generation  of  a  preliminary  version  of  fully  self-consistent  CVB-DFT.  The 
procedure  we  employed  to  accomplish  each  of  these  tasks  is  described  in  this  report,: 
along  with  MCDFT  results  for  computation  of  chemical  properties.  Our  Phase  I  results 
are  highly  encouraging,  demonstrating  the  MCDFT  approaches  are  camable  of  signific¬ 
antly  improved  accuracy  as  compared  to  current  DFT  functionals.  The  N3  or  better 
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pseudospectral  (PS)  algorithm  ensures  that  MCDFT  is  practical  for  reasonable-sized 
systems  as  well.  Although  considerable  additional  work  will  be  required  in  Phase  II 
to  implement  fully  self-consistent  MCDFT  methods,  fund  the  optimum  functional  (or 
functionals),  improve  computational  efficiency,  and  define  protocols  for  utilizing 
GVB  wavefunctions  in  a  localized  region  of  the  molecule,  the  preliminary  software 
clearly  demonstrates  that  the  basic  multiconf igurational  density  functional  theory 
is  sound. 
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The  multiconfigurational  density  functional  theory  (MCDFT)  software  Schrodinger. 
Inc.  has  developed  in  Phase  I  of  this  project  combines  many  of  the  advantages  of  generalized 
valence  bond  (GVB)  and  restricted  configuration  interaction  (RCI)  techniques  with  those  of 
DFT.  We  have  explored  several  avenues  for  improving  the  implementation  of  MCDFT 
methods  and  have  incorporated  the  MCDFT  code  into  a  development  version  of  our  commer¬ 
cially  successful  electronic  structure  program,  PS -GVB  [1], 

The  work  we  have  performed  in  Pha,se  I  has  centered  upon  these  five  tasks:  parti¬ 
tioning  of  the  GVB  two-electron  energy  into  Coulomb,  exchange,  and  intra-pair  terms;  coding 
and  developing  GVB-RCI-DFT  (post-SCF  DFT  on  a  GVB-RCI  density);  inclusion  of  self¬ 
interaction-corrected  (SIC)  density  functionals;  optimization  of  hybrid  MCDFT  methods;  and 
generation  of  a  preliminary  version  of  fully  self-consistent  GVB-DFT.  The  procedure  we 
employed  to  accomplish  each  of  these  tasks  is  described  in  this  report,  along  with  MCDFT 
results  for  computation  of  chemical  properties. 

Our  Phase  I  results  are  highly  encouraging,  demonstrating  that  MCDFT  approaches 
are  capable  of  significantly  improved  accuracy  as  compared  to  current  DFT  functionals.  The 
or  better  scaling  with  basis  set  size  for  PS-GVB’s  methods  resulting  from  its  use  of  the 
pseudospectral  (PS)  algorithm  ensures  that  MCDFT  is  practical  for  reasonable-sized  systems 
as  well.  Although  considerable  additional  work  will  be  required  in  Phase  II  to  implement  fully 
self-consistent  MCDFT  methods,  find  the  optimum  functional  (or  functionals),  improve 
computational  efficiency,  and  define  protocols  for  utilizing  GVB  wavefunctions  in  a  localized 
region  of  the  molecule,  the  preliminary  software  clearly  demonstrates  that  the  basic  multicon¬ 
figurational  density  functional  theory  is  sound. 

Accomplishments  and  New  Findings 

A.  Summary  of  MCDFT  Methods 

The  theoretical  approach  for  MCDFT  methods  originates  with  a  GVB  or  GVB-RCI 
wavefunction,  which  adequately  represents  the  static  correlation  effects  necessary  to  obtain 
the  correct  features  of  a  potential  surface,  such  as  left-right  correlation  and  proper  spin 
coupling.  In  the  GVB  approach,  each  bond  or  other  electron  pair  is  described  by  two  non- 
orthogonal  orbitals,  whose  contributions  to  the  bond  description  are  obtained  variationally. 
Because  PS-GVB  has  a  high-quality  automated  initial  guess  for  the  wavefunction  [2]  and  fast, 
reliable  convergence  algorithms  [3],  its  GVB  module  is  both  highly  efficient  and  straightfor¬ 
ward  to  use.  The  GVB-RCI  program  within  PS-GVB  generates  a  correlated  wavefunction 
from  intra-pair  excitations  of  the  GVB  reference  wavefunction,  using  a  highly  effective 
contraction  procedure  to  reduce  the  length  of  the  Cl  expansions  [4].  The  program  employs  the 
pseudospectral  method  to  speed  up  integral  evaluation  by  reducing  the  scaling  of  the  evalua¬ 
tion  of  each  Coulomb  or  exchange  operator  in  basis  function  space  from  to  N-\  and  solving 
for  the  necessary  matrix  elements  with  a  fast  two-index  transform  rather  than  the  expensive 
four-index  transform  required  in  traditional  ab  initio  codes.  The  program  also  systematically 
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includes  only  the  most  important  configurations  to  make  the  calculation  more  practical,  with 
minimal  loss  of  accuracy  relative  to  the  fully  uncontracted  expansion.  The  internal  contraction 
scheme  used  restricts  the  number  of  Cl  coefficients  in  the  RCI  calculation  to  ~n^,  where  n  is 
the  number  of  GVB  pairs,  yet  is  in  excellent  agreement  with  a  fully  uncontracted  Cl  which 
would  contain  2"n^  Cl  coefficients  (the  number  of  uncontracted  determinants). 

To  add  dynamic  correlation  energy,  E^,  one  can  simply  use  the  GVB  or  GVB-RCl 
charge  density  in  DFT  correlation  functionals  or,  as  a  more  sophisticated  approach,  in  hybrid 
DFT  methods,  which  include  both  GVB  or  GVB-RCl  “exact”  exchange  and  density  func¬ 
tionals  representing  exchange  as  well  as  correlation  functionals.  As  a  modification  to  Becke’s 
theory  [5]  of  hybrid  density  functional  methods,  we  have  also  experimented  with  scaling  an 
intra-pair  GVB  term  (described  below)  as  well  as  scaling  the  DFT  exchange  energy 
and  the  non-local  contribution  to  the  correlation  energy  During  Phase  I,  we  have 

generated  some  preliminary  results  for  the  hybrid  parameters  and  the  parametrized  energy 
partitioning  required  for  these  methods.  We  have  also  begun  to  implement  a  self-consistent 
version  of  GVB-DFT,  as  detailed  in  Section  F. 

Multiconfigurational  density  functional  theory  leads  to  improved  accuracy  at  quite 
reasonable  computational  costs,  and  it  can  be  expected  to  be  partieularly  useful  for  various 
problems  that  are  not  easily  studied  with  standard  ab  initio  and  density  functional  methods. 
For  instance,  transition  metal  chemistry  is  known  to  have  large  multireference  effects,  which 
may  explain  the  marginal  performance  of  DFT  on  these  systems  [6].  The  use  of  an  MCSCF 
reference  also  allows  one  to  treat  a  number  of  electronic  systems  in  which  more  than  one 
determinantal  wavefunction  is  required  even  in  zeroth  order.  For  example,  symmetry  and  spin 
eigenstates  of  many  of  the  states  of  a  system  as  simple  as  O2  require  a  multideterminental 
expansion  which  cannot  be  obtained  within  the  context  of  traditional  density  functional 
theory.  Transition  states  of  so-called  “symmetry  forbidden”  reactions  or  the  calculation  of 
potential  curves  near  avoided  crossings  involves  mixing  of  several  determinantal  functions 
which  again  cannot  be  represented  by  density  functional  theory  alone. 

B.  Partitioning  of  the  GVB  Two-Electron  Energy  into  Coulomb, 
Exchange,  and  Intra-Pair  Terms 

The  key  to  the  performance  of  our  MCSCF-DFT  methods  largely  lies  in  the  ability  of 
the  ab  initio  method  to  robustly  describe  dominant  static  correlation  effects  and  the  excellent 
scaling  of  the  combined  methods  with  basis  set  size.  Therefore,  as  part  of  our  implementation 
of  self-consistent  MCDFT  methods,  we  have  performed  some  of  the  work  necessary  to 
improve  the  treatment  of  the  GVB  reference  wavefunction. 

The  partitioning  of  the  two-electron  energy  into  Coulomb  and  exchange  components 
is  not  so  uniquely  defined  for  GVB  or  GVB-RCl  wave  functions  as  it  is  for  Hartree-Fock,  and 
in  course  of  our  Phase  I  work  we  considered  different  possible  schemes.  The  issue  is  moot,  of 
course,  for  purely  “exact-exchange”  methods,  where  the  DFT  correlation  energy  for  the  final 
GVB-RCl  density  is  simply  added  to  the  complete  final  GVB-RCl  energy.  However,  to  design 
an  accurate  hybrid  GVB-DFT  method,  in  which  exchange  and  correlation  contributions  from 
both  GVB  and  DFT  are  scaled  and  combined,  it  is  quite  important  that  physically  related 
terms  are  first  identified  in  the  GVB  and  DFT  energy  expressions. 
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Since  the  classical  Coulomb  energy  of  the  electron  density  can  be  unambiguously 
defined  in  both  the  GVB  and  DFT  approaches,  it  is  natural  to  separate  this  from  the  remaining 
(nonclassical)  exchange-correlation  contributions  to  the  two-electron  energy;  this  is  the 
normal  approach  in  DFT  theory,  but  not  in  MCSCF  methods  like  GVB.  The  GVB  energy  can 
be  rewritten  so  that  interactions  between  the  orbitals  of  a  given  GVB  pair  and  all  other  orbitals 
are  described  by  “mean-field”  Coulomb  and  exchange  operators,  and  The  mean-field 
Coulomb  operator  thus  defined  corresponds  to  the  Coulomb  operator  of  DFT  theory.  The 
expectation  values  for  these  operators  are  the  classical  Coulomb  and  mean-field  exchange 
energies 


OCC  1 


p(^i)p(/S) 


and 
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i.i  =  1 


/»./■ 


=  - 1  m, 


ij=  I 


where  i  and  /  run  over  all  occupied  orbitals  and  /j-  is  the  occupation  number  of  orbital  i.  For 
closed-shell  systems,  the  remaining  two-electron  energy  consists  entirely  of  terms  local  to 
each  GVB  pair,  and  the  sum  of  these  is  the  intra-pair  “exchange”  energy, 


=  I  cj,  ( 1  -  c,;„)  ,,,,  +  cf,„  ( I 

P 


where  and  C^/,  are  the  Cl  coefficients  for  the  natural  orbitals  and  of  GVB  pairp.  It 
is  important  to  remember  that  also  describes  the  correlation  energy  of  the  GVB  wave- 
function. 

With  these  definitions,  the  complete  energy  expression  for  the  general  open-shell 
(high-spin)  GVB  wavefunction  is 


'GVB 


-  ^nuc  +  ^ f-vb  +  ^open 


where  contains  remaining  non-mean-field  exchange  interactions  among  the  open-shell 
orbitals.  The  complete  exchange-correlation  energy  of  the  GVB  wavefunction  is  therefore  K  ^ 
+  purpose  of  comparing  with  the  DFT  energy  expression, 

corresponds  roughly  to  the  Hartree-Fock  exchange  energy,  while  contains  components  of 
both  exchange  and  correlation. 

To  reiterate,  the  point  of  separating  the  non-Coulomb  part  of  the  two-electroti  energy 
in  this  manner  is  to  try  to  isolate  the  exchange-like  and  correlation-like  components  of  the 
GVB  energy,  so  that  effective  hybrid  methods  can  be  defined  that  minimize  the  double¬ 
counting  of  exchange  and  correlation  when  DFT  functionals  are  introduced.  The  aforemen¬ 
tioned  scheme  represents  just  one  possible  separation.  In  fact,  the  parameterization  results 
obtained  in  Phase  I,  while  quite  good,  often  demanded  that  the  and  components  be 
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scaled  almost  equally;  this  indicates  to  us  that  the  and  defined  here  may  not  be  phys¬ 
ically  independent  and  that  a  different  separation  may  provide  more  flexibility.  For  instance, 
there  are  still  intra-pair  exchange  terms  in  can  be  moved  into  We  will 

examine  the  effect  of  this  and  other  partitionings  in  Phase  II. 

C.  Post-SCF  DFT  on  a  GVB-RCI  Density 

During  Phase  I,  we  implemented  GVB-RCI-DFT  within  PS-GVB,  an  effort  that 
required  computing  RCI  densities  for  a  post-SCF  DFT  analysis.  Preliminary  results  for  post- 
SCF  correlation-only  DFT  treatments  showed  a  fairly  low  average  error  of  5.9  kcal/mol  in 
predictions  of  experimental  atomization  energies  of  closed-shell  systems. 

We  have  also  used  the  GVB-RCI-DFT  software  to  evaluate  conformational  energies 
calculated  by  performing  a  post-SCF  DFT  analysis  on  the  RCI  density  using  various  local 
correlation  functionals.  Table  1  lists  the  results  we  have  obtained  thus  far  and  compares  them 
with  values  obtained  using  GVB-RCI  alone.  In  most  cases  the  DFT  correction  results  is  an 
improvement  of  the  RCI  reference  value.  Given  that  we  have  not  parametrized  the  RCI 
exchange  and  had  not  yet  included  self-interaction  corrections,  these  results  suggest  that  the 
RCI-DFT  method  is  sound. 

In  Table  2,  we  present  timings  for  GVB-RCI-DFT  calculations  on  systems  with  up  to 
287  basis  functions  and  21  GVB  pairs,  using  a  cc-pVTZ(-f)  basis  set.  To  expedite  the  GVB 
calculations,  we  first  converged  6-3 IG**  wavefunctions,  then  used  those  as  initial  guesses  for 
the  cc-pVTZ  GVB  runs.  The  GVB  calculation  is  clearly  the  rate-limiting  factor  of  these  calcu¬ 
lations;  this  effect  will  be  mitigated  by  our  strategies  to  optimize  GVB  performance  in  Phase 
II. 


The  overall  .scaling  of  the  GVB-RCI-DFT  method  is  ,  which  is  far  superior  to  any 
other  common  MCSCF  procedure  such  as  CASSCF.  Given  these  timings,  routine  application 
of  GVB-RCI-DFT  to  systems  with  on  the  order  of  500  basis  functions  and  40  GVB  pairs  (80 
correlated  electrons)  is  feasible  on  workstations  with  reasonable  throughput.  Furthermore,  not 
all  GVB  pairs  need  to  be  included  in  a  given  problem,  an  advantage  which  further  extends  the 
range  of  applicability  of  these  methods. 

We  are  continuing  to  address  the  issue  of  partitioning  the  GVB-RCI  exchange  term 
into  mean  field  and  intra-pair  contributions  in  order  to  extract  the  analogue  of  the  intra-pair 
GVB  term  described  in  Section  B.  Although  the  RCI  energy  expression  is  considerably  more 
complicated  than  the  GVB  expression,  the  energy  can  be  partitioned  into  Coulomb  and 
exchange  terms  as  well  as  intra-  and  inter-pair  terms.  We  anticipate  that  for  best  results  with 
hybrid  methods,  we  will  then  need  to  scale  this  intra-pair  RCI  term,  as  we  did  for  GVB -DFT. 

D.  Correlation  Functional  Modifications 

We  have  found,  as  have  others  [7,  8,  9],  that  a  correlation  correction  calculated  with 
standard  DFT  correlation  functionals  using  the  GVB  density  leads  to  an  over-correction  of  the 
correlation  energy.  For  example,  an  atomization  calculation  of  C2H4  has  a  GVB  value  of 
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Molecule 

RCI 

RCI-VWN 

RCI-VWN5 

RCI-PZ81 

Exp. 

cyclohexanol 

1.04 

0.93 

0.49 

0.47 

0.52 

piperidine 

0.86 

0.87 

0.87 

0.87 

0.40 

2,3-dimethylbutane 

0.03 

0.04 

0.03 

0.03 

0.17 

N-methylform  amide 

0.97 

1.10 

1.09 

1.08 

1.45 

butane 

1.16 

0.97 

1.00 

1.00 

0.75 

methyl  acetate 

8.6! 

8.71 

8.73 

8.71 

7.5-8.5 

methyl  ethyl  ether 

1.71 

1.55 

1.58 

1.58 

1.50 

acrolein 

1.82 

1.73 

1.75 

1.75 

2.00 

isopropylamine 

0.31 

0.43 

0.41 

0.41 

0.45 

propionaldehyde 

0.64 

0.83 

0.80 

0.80 

0.95 

methyl  formate 

4.97 

5.36 

5.31 

5.30 

4.75 

ethyl  ether 

1.87 

1.53 

1.60 

1.60 

1.10 

1-butene 

0.72 

0.41 

0.45 

0.45 

0.53 

2-butene 

1.89 

1.68 

1.72 

1.71 

1.00 

butanone 

1.09 

1.27 

1.24 

1.24 

1,15 

1,3-butadiene 

2.49 

2.45 

2.47 

2.46 

2.49 

methyl  vinyl  ether 

0.38 

0.81 

0.75 

0.75 

1.15 

cyclohexamine 

1.42 

0.76 

0.89 

0.89 

1.10 

Table  1:  Conformational  energies  in  kcal/mol,  evaluated  with  GVB-RCI  only  (RCI)  and  by 
performing  a  post-SCF  DFT  analysis  of  the  GVB-RCI  self-consistent  wavefunction  using 
various  local  correlation  functionals. 


Molecule 

^bas 

N  . 

pair 

6-31G** 

cc-pVTZ 

RCI  Int. 

RCI  Energy 

DFT 

Total 

Methylcyclohexane 

287 

21 

76 

550 

110 

6 

142 

884 

Cyclohexane 

246 

18 

43 

340 

68 

3 

86 

540 

Methylvinylether 

146 

12 

12 

83 

10 

0.5 

13 

119 

Table  2:  CPU  times  (minutes)  for  pseudospectral  GVB  initial  wavefunction  calculations  with  the 
6-31G**  basis  and  with  the  cc-pVTZ(-f)  basis  using  the  6-31G*'^  results  as  initial  guesses,  RCI 
integral  generation,  RCI  energy  solver,  and  DFT  correlation  energy  treatment  of  the  GVB-RCI 
density.  All  calculations  performed  on  a  single  IBM-SP2  390  thin  node. 
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461.4  kcal/mol  and  a  VWN  correlation  functional  corrected  value  of  581.1  kcal/mol,  in 
comparison  to  the  correct  atomization  energy  of  531.9  kcal/  mol. 

There  are  two  established  methods  for  dealing  with  this  over-correlation  that  need  to 
be  implemented  and  tested  within  the  GVB-DFT  and  GVB-RCTDFT  formalisms:  modifying 
the  spin  part  of  the  density  interactions,  and  modifying  the  total  spin-space  density  that  enter 
into  the  correlation  functionals.  Although  we  originally  slated  work  on  both  of  these  methods 
for  Phase  II,  we  have  partially  completed  programming  and  testing  one  of  these  approaches, 
the  inclusion  of  “self-interaction-corrected”  (SIC)  functionals  [7,  10]  within  the  code.  These 
functionals  assume  that  the  dominant  correlation  effect  is  between  electrons  of  opposite  spins 
and  thus  subtract  out  the  correlation  from  like  spins,  writing  the  DFT  correlation  correction  to 
the  energy  as: 

E,  =  prp(r)ejp^^  (;-),pp(r)]  -  jr/rp„  (r)  £,  lp„  (r) ,  0]  -prpp  (r)  £,.  [0,  pp  (r)  ] 

where  and  pp  denote  the  aa  and  (3(3  parts  of  the  total  spin  density  p  =  p^^  -h  pp  and 
e^[poj(r),  pp(r)]  is  a  standard  correlation  functional.  This  approach  could  be  further  parame¬ 
trized  by  scaling  (with  a  factor  of  0  to  1)  the  part  of  the  self-interaction  which  is  subtracted  in 
the  above  equation.  Note  that  without  the  SIC  correction,  the  H  atom  has  an  unphysical  non¬ 
zero  DFT  correlation  energy.  Part  of  the  rationale  for  this  approach  is  that  the  treatment  of  ab 
initio  (HF  or  GVB)  exchange  interactions  is  providing  for  a  large  fraction  of  the  aa,  PP  corre¬ 
lations  and  hence  this  correlation  should  be  removed  from  as  above. 

Some  preliminary  results  for  atomization  energy  calculations  generated  using  the  new 
SIC  code  are  included  in  Section  E.2  below.  In  Phase  II,  we  will  continue  to  study  the  perfor¬ 
mance  of  such  SIC  modifications  for  various  correlation  functionals  and  in  conjunction  with 
the  exchange  scaling  schemes  discussed  in  Section  E. 

E.  Optimization  of  Hybrid  MCDFT  Methods 

E.l.  Conformational  Energy  Calculations 

The  “three  parameter”  method  that  was  developed  by  Becke  [5]  and  is  now  widely 
used  for  density  functional  theory  calculations  was  designed  to  give  good  self-consistent  DFT 
results  by  combining  optimal  amounts  of  “exact”  Hartree-Fock  and  density  functional  terms. 
The  parameters  were  chosen  by  performing  least-squares  fitting  to  obtain  the  values  that 
yielded  the  best  results  for  chemical  properties  of  several  dozen  molecules.  Popular  hybrid 
methods  that  use  these  thi'ee  parameters  include  B3LYP. 

Because  Becke’s  parametrization  was  designed  to  allow  the  DFT  treatment  to  account 
for  effects  neglected  by  Hartree-Fock,  we  did  not  expect  the  parameters  to  be  ideal  for  post- 
SCF  DFT  analysis  of  GVB  wavefunctions  or  for  self-consistent  GVB-DFT.  In  Phase  I,  we 
performed  a  preliminary  analysis  of  what  these  parameters  should  be  to  give  best  results  for 
GVB-DFT  and  found  that  re-fitting  the  parameters  gave  a  dramatic  improvement  for  post- 
SCF  DFT  calculations  on  GVB  wavefunctions,  regardless  of  the  local  and  non-local  correla¬ 
tion  functionals  used. 
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Our  preliminary  GVB-DFT  parametrization  was  designed  to  minimize  RMS  errors  for 
a  set  of  17  conformational  energy  differences,  where  the  cases  and  their  experimental  values 
were  drawn  from  among  papers  concerning  calculation  of  relative  conformational  energies  by 
Murphy  et  al.  [11]  and  St.-Amant  et  al.  [12].  The  conformational  energy  calculations  were 
performed  with  our  electronic  structure  software,  PS-GVB,  using  a  cc-pVTZ(-f)  basis  set,  an 
energy  convergence  criterion  of  l.Ox  10'^  Hartrees,  and  PS-GVB’s  highest-accuracy  setting. 
The  DFT  exchange-correlation  energy  was  calculated  non-self-consistently  using  the  self- 
consistent  GVB  density;  all  GVB  calculations  included  all  possible  sigma  and  pi  bond  pairs. 

The  RMS  error  obtained  for  pure  GVB  conformational  energy  differences,  0.43  kcal/ 
mol,  was  considerably  lower  than  the  HF  or  B3LYP  RMS  errors  of  0.52  and  0.54  kcal/mol, 
and  including  a  treatment  of  the  GVB  densities  with  correlation  functionals  yields  markedly 
improved  RMS  errors  (0.32  kcal/mol  for  GVB-LYP,  for  instance,  an  improvement  of  more 
than  30%  over  the  pure  HF  or  B3LYP  results).  However,  applying  the  hybrid  DFT  method 
B3LYP  to  the  GVB  densities  gave  a  quite  poor  RMS  error  of  0.56  kcal/mol.  Clearly  the 
weights  assigned  to  various  terms  by  Becke’s  three  parameters  were  unsuitable  for  our  GVB- 
DFT  conformational  energy  calculations. 

Results  of  our  Phase  I  parametrization  of  GVB-DFT  hybrid  methods  are  summarized 
in  Table  3.  We  first  performed  a  completely  unrestricted  parameter  fit  to  minimize  the  RMS 
error,  analyzing  coefficients  for  every  GVB  and  DFT  term,  with  no  restrictions  on  the  range  of 
the  coefficients.  During  this  fit,  we  allowed  the  coefficients  for  the  exact  exchange  terms 
and  (described  in  Section  B)  to  vary  separately.  As  we  expected  and  as  Table  4  shows, 
the  resultant  coefficients  were  totally  unphysical,  but  this  fitting  gave  us  an  idea  of  the  best 
result  we  could  possibly  expect  for  an  RMS  conformational  energy  difference  error  for  the 
data  set  in  question,  using  only  GVB  and  post-SCF  DFT  terms:  an  RMS  error  of  0.23  kcal/ 
mol,  considerably  lower  than  the  GVB  and  RMS  error  of  0.43  and  less  than  half  the  B3LYP 
error  of  0.56.  The  terms  from  the  local  correlation  functional  using  the  GVB  density  are 
clearly  the  most  important  for  a  good  fit  to  experimental  data.  The  unrestricted  fit  also  makes 
it  clear  that  the  Slater  local  exchange  functional  [13],  whose  coefficient  is  very  small,  is 
contributing  nothing  useful  to  the  fit:  the  GVB  exact  exchange  dominated  the  calculation  of 
exchange  energy.  Therefore,  for  the  remaining  fittings,  we  eliminated  the  Slater  contribution. 


Pure  SCF 

Post-SCF  DFT  on  GVB  density 

Single  functional 

Parametrized  combination  of  functionals 

GVB 

B3LYP 

VWN 

LYP 

GGA-IIC 

B3LYP 

unrestricted 

VVVN/P86 

PW91/P86 

VWN/GGAIIc 

PW91/GGAIIC 

VWN/LYP 

PW91/LYP 

Error: 

0.43 

0.54 

0.33 

0.32 

0.34 

0.56 

0.23 

0.30 

0.31 

0.31 

0.31 

0.31 

0.30 

0.31 

Table  3:  RMS  errors  in  kcal/mol  for  calculations  of  17  conformational  energy  differences  using 
various  SCF  and  post-SCF  DFT  methods.  Coefficients  used  for  parametrized  combinations  of 
functionals  are  shown  in  Table  4  and  Table  5. 


8 


Fitting 

Conditions 

Slater 

Becke88 

VWN 

PW-91 

LYP 

Perde\v86 

GGA-lIc 

unrestricted 

1.3014 

0.9530 

-0.2031 

-7.0147 

104.19 

-116.42 

-2.6362 

-24.585 

12.282 

1.0000 

0.9477 

0.0000 

0.5497 

1 .0000 

0.7813 

0.5628 

0.5265 

0.0000 

Table  4:  Coefficients  for  various  terms  for  parameter  fitting  to  reduce  RMS  errors  for 
conformational  energy  differences  described  in  the  text.  The  1  fit  also  restricted  all  other 
coefficients  to  be  between  0  and  1. 

When  we  constrained  the  coefficient  to  unity  and  restricted  the  other  coefficients 
to  be  between  0  and  1,  while  still  allowing  all  possible  functionals  to  contribute,  the  RMS 
error  was  0.30,  still  much  lower  than  GVB  or  B3LYP  alone  gave.  However,  as  Table  4  shows, 
the  total  coefficients  for  each  type  of  exchange  or  correlation  term  were  too  high.  In  addition, 
including  a  contribution  from  every  functional  available  would  be  difficult  to  justify  as  a 
general  method. 

Our  next  goal  was  to  restrict  the  parameters  to  reasonable  values  and  definitions 
without  significantly  reducing  the  accuracy  of  the  GVB-DFT  results.  We  followed  the  model 
provided  by  Becke’s  three  parameters,  restricting  the  sum  of  the  parameters  for  total  exact 
exchange  and  local  exchange  to  1  (in  this  case  by  including  all  of  the  GVB  exchange 
energy  and  leaving  out  the  Slater  term),  and  ensuring  that  the  parameters  for  local  correlation 
also  summed  to  1.  Parameters  for  Becke’s  1988  gradient  correction  to  the  exchange  [14]  and 
for  a  non-local  correlation  functional  were  allowed  to  vary  between  0  and  1.  (As  usual,  when 
LYP  was  used  as  a  non-local  correlation  functional,  the  local  correlation  functional’s  contribu¬ 
tion  was  reduced  correspondingly  to  make  up  for  the  local  correlation  contribution  of  the  LYP 
functional.)  We  also  introduced  another  parameter  for  the  GVB  intra-pair  term  which  is 
described  in  detail  in  Section  B.  In  summary,  we  were  fitting  for  the  a,  (3,  and  x  that  would 
give  the  lowest  RMS  energy  error  for  energies  of  the  form: 


„GVB-DFT  „GVB  „GVB  „GVB  n^NLOA,  ,  x  ry  .  ,, 

^  =  ^Coulomb  +  ^xmf  +  + P^.v  [P  J 

+  f  p  (v|/^„,/,)  ]  +  [p  ( V,„./P .  Vp  (V,,,/,)  J 


while  restricting  all  three  parameters  to  be  between  0  and  1 .  We  optimized  these  parameters 
for  six  different  cases,  where  for  each  individual  case,  we  chose  either  VWN  [15]  or  PW-91 
[16]  as  the  local  correlation  functional  and  either  LYP  [17],  Perdew  86  [18],  or  GGA-IIc  [16] 
as  the  non-local  correlation  functional. 

Ideally,  we  wanted  our  results  to  be  nearly  as  accurate  as  the  0.30  RMS  error  we  had 
obtained  with  the  minimal  constraints  described  above,  and  to  remain  relatively  stable  with 
our  choice  of  functional,  particularly  in  regard  to  the  local  correlation  functional,  which  is  not 
scaled  by  any  variable  parameter.  We  also  felt  the  P  parameter  should  be  low,  since  the  Becke 
gradient  correction  is  really  meant  to  address  the  limitations  of  the  Slater  exchange  functional 
and  should  therefore  be  unnecessary  when  only  the  accurate  GVB  exchange  is  used  instead. 
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Finally,  we  wanted  the  GVB  intra-pair  scaling  parameter  to  be  consistent  for  the  various  func¬ 
tional  choices. 

As  Table  5  indicates,  all  of  these  conditions  were  met:  the  RMS  error  ranged  from  0.30 
to  0.3 1  kcal/mol  over  all  six  cases,  the  value  obtained  for  (3  in  every  case  was  0,  and  the  value 
for  X  was  always  0.94.  This  preliminary  work  should  provide  a  solid  foundation  for  further 


VWN/ 

P86 

PW91/ 

P86 

VWN/ 

GGAIIc 

PW91/ 

GGAIIc 

VWN/ 

LYP 

PW91/ 

LYP 

a 

0.94 

0.94 

0.94 

0.94 

0.94 

0.94 

p 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

X 

0.25 

0.33 

0.29 

0.40 

0.42 

0.54 

RMS  Error 

0.31 

0.31 

0.31 

0.30 

0.31 

0.30 

MA  Error 

0.24 

0.24 

0.24 

0.23 

; _ i 

0.23 

0.23 

Table  5:  Parameters  fitted  for  various  combinations  of  functionals  for  post-SCF  DFT 
calculations,  and  the  resultant  RMS  and  mean  absolute  errors  (in  kcal/mol)  for  conformational 
energy  difference  calculations.  Parameters  are  described  in  the  text. 

parametrization  efforts  in  Phase  II. 

E.2.  Atomization  Energy  Calculations 

Our  Phase  I  calculations  of  GVB-DFT  atomization  energies  for  44  molecules  from  the 
G2  data  set  [19]  further  demonstrate  the  promise  of  MCDFT  methods,  although  much  work 
still  remains  to  be  done  in  Phase  II  to  optimize  hybrid  method  parameters  and  other  aspects  of 
GVB-DFT  and  GVB-RCI-DFT  calculations.  All  of  our  preliminary  atomization  energy  calcu¬ 
lations  described  here  used  the  cc-pVTZ(-f)  basis  set  and  PS-GVB’s  highest  accuracy  setting, 
starting  with  the  standard  initial-guess  GVB  wavefunction.  For  the  first  set  of  calculations 
described  here,  those  without  self-interaction-corrected  functionals,  GVB  pairs  were  defined 
for  all  molecular  bonds,  but  not  lone  pairs,  while  atoms  were  treated  at  the  ROHF  level;  for 
the  calculations  with  self-interaction  corrections,  GVB  lone  pairs  were  also  included. 

We  first  evaluated  the  effects  of  various  functionals  in  the  post-SCF  DFT  treatment  of 
HF  and  GVB  wavefunctions.  The  use  of  local  and  non-local  DFT  correlation  functionals 
uniformly  gave  a  dramatic  improvement  over  the  straight  HF  and  GVB  atomization  energies, 
as  is  shown  in  Table  6.  However,  the  unsealed  inclusion  of  the  Slater  local  exchange  func¬ 
tional  and  Becke  non-local  gradient  correction  to  the  exchange  (BLYP,  for  example)  actually 
made  the  GVB-based  results  worse,  although  the  HF-based  atomization  energies  again 
improved  substantially.  We  did  ultimately  obtain  higher  quality  results  for  the  GVB-based 
treatment  when  including  DFT  exchange  functionals,  but  only  by  retaining  some  fraction  of 
the  GVB  exchange  energy  in  hybrid  methods  or  by  using  SIC  functionals.  Some  of  our  results 
for  hybrid  methods  are  included  in  Table  6  for  comparison;  the  methods  are  explained  in 
detail  below. 
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SCF 

wavefunction 

DFT 

correlation 

functional(s) 

DFT 

exchange 

functiona](s) 

DFT 

hybrid 

method 

atomization 

energy 

error 

HF 

— 

— 

— 

80.3 

GVB 

— 

— 

— 

49.7 

HF 

VWN  (local)  only 

— 

— 

26.0 

GVB 

VWN  (local)  only 

— 

— 

23.7 

HF 

PW91  (local)  only 

— 

— 

28.6 

GVB 

PW91  (local)  only 

— 

— 

22.8 

HF 

GGA-II  (non-local)  only 

— 

— 

25.2 

G 

Co 

HF 

LYP  (non-local)  only 

— 

““ 

25.6 

Co 

cS 

GVB 

GGA-II  (non-local)  only 

— 

— 

17.2 

GVB 

LYP  (non-local)  only 

— 

— 

18.4 

i 

^  s- 

J  § 

HF 

GGA-II  (non-local)  only 

GGA-II 

— 

5.1 

HF 

LYP  (non-local)  only 

Slater/Becke 

— 

5.3 

HF 

VWN  (local)  /  LYP  (non-local) 

Slater/Becke 

B3LYP 

5.0 

k, 

GVB 

GGA-II  (non-local)  only 

GGA-II 

— 

26.4 

GVB 

LYP  (non-local)  only 

Slater/Becke 

— 

30.7 

GVB 

VWN  (local)  /  LYP  (non-local) 

Slater/Becke 

B3LYP 

19.8 

GVB 

VWN  (local)  /  Perdew86  (non-local) 

Slater/Becke 

new 

4.06 

G  c 

Co  -Ci 

g  ^ 

GVB 

PW91  (local)  /  Perdew  86  (non-local) 

Slater/Becke 

new 

4.11 

GVB 

VWN  (local)  /  GGA-IIc  (non-local) 

Slater/Becke 

new 

3.36 

GVB 

PW91  (local)  /  GGA-IIc  (non-local) 

Slater/Becke 

new 

3.73 

"O 

GVB 

VWN  (local)  /  LYP  (non-local) 

Slater/Becke 

new 

3.45 

GVB 

PW9I  (local)  /  LYP  (non-local) 

Slater/Becke 

new 

3.36 

Table  6:  Full  comparison  of  mean  absolute  errors  (in  kcal/mol)  for  atomization  energy 
calculations  on  50  molecules  from  the  G2  data  set  using  various  ab  initio  methods.  All  DFT 
calculations  were  post-SCF  treatments  of  the  DFT  energy  using  the  density  from  the  given  self- 
consistent  wavefunction.  New  hybrid  methods  are  explained  later  in  this  section.  Entries  of  — 
mean  no,  none,  or  not  relevant. 


Hybrid  methods  we  considered  contained  various  proportions  of  the  GVB  mean-field 
exchange  energy,  the  GVB  intra-pair  energy,  the  Slater  local  exchange  functional,  the  Becke 
non-local  exchange  functional,  the  VWN  and  Perdew-Wang  1991  local  correlation  func¬ 
tionals,  and  the  Perdew  1986,  GGA-II  and  LYP  non-local  correlation  functionals.  We  first 
tried  to  optimize  only  the  a,  |3,  and  %  parameters  described  in  Section  E.l,  leaving  out  the 


Slater  exchange  functional  term  (and  including  all  of  the  exact  exchange),  and 

parametrizing  for  various  combinations  of  one  local  and  one  non-local  correlation  functional. 
The  resultant  RMS  errors  were  reasonable,  but  somewhat  disappointing,  ranging  from  15-20 
kcal/mol  depending  on  the  functional  choices. 

Noting  that  the  large  RMS  error  of  53.4  kcal/mol  for  the  GVB-only  calculations  might 
indicate  problems  with  the  GVB  reference  wavefunctions,  we  next  allowed  a  non-zero  contri¬ 
bution  for  the  Slater  exchange,  subtracting  off  a  corresponding  amount  of  the  exact  exchange. 
The  intra-pair  term  was  still  allowed  to  vary  independently.  The  total  GVB-DFT  energy  was 
thus  taken  to  be 


„GVB-DFT 

E 


GVB  GVB  .  ^Sillier  GVB 

~  CouUnnb'^  ^  X  kiiif  ^  ^.v  P  ^^iiitra  -  pair 

+  [p  ,  V p  „)  J  +  [p  ] 

+  [ P  ( VI/,,,/,) ,  Vp  ( y ] 


where  we  constrained  all  four  parameters  (c,,  a,  P,  and  x)  to  be  between  zero  and  one.  The 
RMS  errors  improved  dramatically  for  all  hybridization  schemes  examined,  as  shown  in 
Table  7. 


Pure  SCF 

Post-SCF  DFT  on  GVB  density 

Single  functional 

Parametrized  combination  of  functionals 

GVB 

B3LYP 

NAVA 

LYP 

GGA-IIC 

B3LYP 

VWN/P86 

PW91/P86 

VWN/GGA-IIc 

PW91/GGA-IIC 

VWN/LYP 

PW91/LYP 

RMS  Error: 

53.4 

3.5 

27.8 

23.7 

22.5 

19.8 

4.95 

5.00 

4.57 

4.66 

4.43 

4.46 

Table  7:  RMS  errors  in  kcal/mol  for  atomization  energy  calculations  for  44  molecules  using 
various  SCF  and  post-SCF  DFT  methods.  Coefficients  used  for  parametrized  combinations  of 
functionals  are  shown  in  Table  8. 


The  optimized  parameters  for  the  various  hybrid  functionals,  shown  in  Table  8,  show  a 
number  of  notable  trends.  The  a  and  c,  coefficients  for  the  GVB  mean-field  exchange  and 
intra-pair  terms  are  quite  similar  for  each  individual  method,  as  seen  also  in  the  conforma¬ 
tional  energy  difference  study;  this  may  indicate  that  our  current  separation  of  these  compo¬ 
nents  has  not  yet  isolated  two  physically  distinct  contributions  to  the  total  GVB  exchange- 
correlation  energy.  Unlike  the  conformational  energy  study,  here  a  large  proportion  of  the 
Slater  exchange  energy  was  required  to  achieve  best  results,  which  previously  were  only 
comparable  to  that  in  the  standard  B3LYP  parameterization.  (We  note  that  Becke’s  three- 
parameter  method  was  optimized  to  produce  good  atomization  energies  using  molecules  from 
this  same  data  set,  which  probably  explains  why  its  performance  was  fairly  bad  for  the  confor¬ 
mational  energies  we  considered  yet  quite  good  for  these  atomization  energies.) 
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B3LYP 

VWN/ 

P86 

PW91/ 

P86 

VWN/ 

GGAIIc 

PW91/ 

GGAIIc 

VWN/ 

LYP 

PW9I/ 

LYP 

Cx 

0.20 

0.16 

0.18 

0.23 

0.22 

0.30 

0.30 

a 

0.20 

0.18 

0.19 

0.24 

0.24 

0.31 

0.36 

p 

0J2 

0.30 

0.20 

0.31 

0.20 

0.24 

0.22 

X 

0.81 

0.00 

0.00 

0.29 

0.20 

0.73 

1 .00 

RMS  Error 

19.8 

4.96 

5.00 

4.57 

4.66 

4.43 

4.46 

MA  Error 

17.5 

4.06 

4.11 

3.36 

3.73 

3.45 

3.36 

Table  8:  Parameters  fitted  for  various  combinations  of  functionals  for  post-SCF  DFT 
calculations,  and  the  resultant  RMS  and  mean  absolute  errors  (in  kcal/mol)  for  atomization 
energy  calculations.  Parameters  are  described  in  the  text. 


With  the  Phase  I  implementation  of  a  self-interaction-corrected  VWN  local  correlation 
functional,  as  described  in  Section  D,  we  reevaluated  atomization  energies  for  a  smaller  set  of 
16  of  the  G2  molecules  that  contained  only  hydrogen,  cai'bon,  nitrogen,  and  oxygen.  For  these 
calculations,  we  also  included  GVB  lone  pairs  as  well  as  GVB  bond  pairs.  Table  9  shows  that 


1 

SCF 

wavefunction 

DFT 

correlation 

functional(s) 

DFT 

exchange 

functional(s) 

DFT 

hybrid 

method 

mean 

absolute 

error 

RMS 

error 

& 

CO 

GVB 

none 

none 

none 

62.9 

66.4 

SCF  B3LYP 

VWN  /LYP 

Slater  /  Becke 

B3LYP 

1.4 

2.0 

q 

GVB 

VWN  (SIC)  only 

none 

none 

14.2 

19.0 

& 

CO 

GVB 

VWN  (SIC) /LYP 

Slater  /  Becke 

B3LYP 

20.2 

21.9 

GVB 

VWN  (SIC)  /  LYP 

Slater  /  Becke 

new  (SIC) 

1.3 

1.8 

Table  9:  Full  comparison  of  mean  absolute  and  RMS  errors  (in  kcal/mol)  for  atomization  energy 
calculations  on  16  molecules  from  the  G2  data  set  using  various  methods,  where  the  VWN  local 
correlation  functional  was  self-interaction-corrected  (SIC)  where  indicated.  All  post-SCF  DFT 
calculations  used  the  self-consistent  GVB  density,  as  indicated.  New  parameters  used  for  last 
calculation  listed  were  c^=0.12,  a=0.25,  (3=0.42,  and  x=0.08,  where  these  parameters  are  as 
described  earlier. 

B3LYP  performs  particularly  well  for  these  molecules,  while  GVB’s  performance  is  not  espe¬ 
cially  good;  nevertheless,  the  highest  accuracy  obtained  is  for  a  post-SCF  DFT  treatment  of 
the  GVB  density,  even  for  this  preliminary,  non-self-consistent  version  of  the  method,  an 
extremely  encouraging  result. 
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We  are  confident  that  our  Phase  II  improvements,  which  will  include  modification  of 
the  GVB  reference  wavefunction,  further  development  of  the  self-interaction-correcting  DFT 
functional  code,  and,  especially,  full  implementation  and  parametrization  of  self-consistent 
GVB-DFT  computations,  will  reduce  errors  still  further.  We  found  that  the  self-consistent  HF- 
B3LYP  mean  absolute  error  for  these  atomization  energy  calculations  is  about  half  that 
resulting  from  a  post-SCF  B3LYP  treatment  of  HF  wavefunctions,  demonstrating  the  impor¬ 
tance  of  a  self-consistent  treatment. 

F.  Self-Consistent  GVB-DFT:  Theory  and  Preliminary  Results 

In  order  to  efficiently  optimize  accurate  molecular  structures  with  the  GVB-DFT  and 
GVB-RCI-DFT  methods,  and  to  calculate  wavefunction-dependent  properties  such  as  multi¬ 
pole  moments,  polarizabilities,  hyperpolarizabilities,  it  is  necessary  that  the  wavefunction  be 
calculated  self-consistently. 

Although  we  originally  planned  to  implement  a  self-consistent  GVB-DFT  method 
only  during  Phase  II,  our  progress  in  Phase  I  encoui’aged  us  to  investigate  the  advantages  of  a 
completely  self-consistent  treatment  as  soon  as  possible.  Therefore,  we  have  completed  a 
preliminary  version  of  this  software,  specialized  to  closed-shell  systems  and  “exact- 
exchange”  methods,  in  which  the  DFT  exchange-correlation  functional  is  used  in  the  refine¬ 
ment  of  the  GVB  orbitals  and  Cl  coefficients.  To  our  knowledge,  this  is  the  only  existing  soft¬ 
ware  with  this  capability. 

In  a  GVB  wavefunction,  each  GVB  electron  pair,  p,  is  described  by  a  pair  of  natural 
orbitals  and  and  their  Cl  coefficients,  and  (where  ^  -r  =  I  ).  To  opti¬ 
mize  the  GVB  wavefunction,  a  separate  Fock  matrix  is  defined  for  the  core  (non-GVB) 
orbitals  and  for  each  GVB  natural  orbital  and  open-shell  orbital  [20].  The  addition  of  the  DFT 
exchange-correlation  energy  to  the  GVB  energy  expression  modifies  the  normal  GVB  Fock 
matrices  by  adding  the  DFT  exchange-correlation  potential  to  each  Fock  matrix;  the  Fock 
matrix,  F-,  for  orbital  i  is 

„GVB-DFT  „GVB 

where  f-  is  the  occupation  of  the  orbital  i  ij]  =  C--  for  GVB  natural  orbital  \|/-)  and  is  the 
DFT  exchange-correlation  potential. 

The  Cl  coefficients  are  optimized  at  each  iteration  to  minimize  the  total  energy  of  the 
wavefunction.  Normally  this  can  be  done  by  solving  a  simple  quadratic  equation  for  the  coef¬ 
ficients  for  each  GVB  pair  in  turn,  and  iterating  over  the  complete  set  of  GVB  pairs  until  self- 
consistency  is  reached.  When  the  complete  GVB  two-electron  energy  is  included  in  the  GVB- 
DFT  energy,  the  equations  for  the  Cl  coefficients  are  only  slightly  modified  in  the  presence  of 
the  DFT  exchange-correlation  terms.  When  the  GVB  exchange  energy  is  scaled,  however,  as 
in  a  hybrid  method  such  as  B3LYP,  a  fourth-order  equation  must  be  solved  for  the  Cl  coeffi¬ 
cients  for  each  pair.  This  equation  and  its  component  teiTns  take  the  form 
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Here,  x  is  a  coefficient  that  scales  the  GVB  exchange  energy  (both  and  similar 

equations  are  obtained  when  the  mean-field  and  intra-pair  energies  are  scaled  independently. 
For  x=l,  the  “exact  exchange”  limit,  the  Z  terms  vanish  and  the  usual  quadratic  equation  for  a 
is  obtained.  For  x=0,  both  and  are  completely  removed  and  the  GVB  wavefunction 
collapses  to  a  single  Slater  determinant,  because  the  DFT  exchange-correlation  energy  does 
not  provide  the  intra-pair  interactions  required  to  stabilize  a  multi-configuration  wavefunc¬ 
tion.  Therefore,  any  hybrid  functional  defined  for  use  in  a  self-consistent  GVB-DFT  theory 
must  retain  some  amount  of  the  GVB  intra-pair  exchange  interaction,  or  it  reduces  to  a  canon¬ 
ical  DFT  theory. 

Solution  of  these  equations  is  straightforward,  and  may  be  simplified  by  using  the  root 
of  the  corresponding  quadratic  equation  as  an  initial  guess.  These  equations  also  make  it  clear 
that  if  the  intra-pair  energy,  completely  scaled  away,  then  the  only  the  lower-energy 

natural  orbital  in  each  pair  is  populated  and  the  GVB  wavefunction  collapses  to  a  single-deter¬ 
minant  HF  wavefunction. 

The  extension  of  these  equations  to  open-shell  systems  is  easily  done.  Open  shell 
systems  are  treated  in  DFT  theory  analogously  to  the  unrestricted  Hartree-Fock  method,  with 
separate  exchange-correlation  potentials  and  for  the  alpha  and  beta  spin-orbitals.  To 
use  these  terms  in  the  context  of  the  explicitly  restricted  GVB  theory,  we  will  employ  the  stan¬ 
dard  relationships  to  recombine  them  into  restricted  core  and  open-shell  potentials, 


closed 


\z'‘  -  1 

open  9  a 


is  then  used  in  the  Fock  matrices  for  the  core  and  GVB  pair  orbitals,  while  is 

used  for  the  open-shell  Fock  matrix. 

We  have  implemented  the  self-consistent  GVB-DFT  theory  described  above  for 
closed-shell  systems  and  “exact  exchange”  methods  and  have  obtained  encouraging  results 
for  the  handful  of  cases  that  we  have  examined.  Convergence  of  the  hybrid  GVB-DFT  wave- 
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function  to  self-consistency  appears  to  be  reliable  using  the  standard  GVB-DIIS  [3]  and 
OCBSE  [20]  convergence  schemes  available  in  PS-GVB,  both  when  starting  from  PS-GVB’s 
default  GVB  initial  guess  and  when  starting  from  a  converged  GVB  wavefunction. 

Table  10  and  Table  1 1  show  preliminary  results  from  three  systems  in  our  conforma¬ 
tional  energy  database.  In  these  systems  there  is  little  additional  lowering  of  the  total  GVB- 
DFT  energy,  when  compared  with  the  non-self-consistent  post-SCF  DFT  treatment,  and  no 
impact  on  the  accuracy  of  the  conformational  energy  differences.  When  the  open-shell  version 


Method 

formic  acid 

1-butene 

isoprene 

(c) 

(t) 

(s) 

(c) 

(g) 

(t) 

GVB 

-188.974709 

-188.967843 

-156.346391 

-156.344988 

-194.236185 

-194.238324 

Post-SCF  DFT 

-191.070371 

-191.063307 

-158.945255 

-158.944344 

-197.346524 

-197.348978 

SCF  GVB-DFT 

-191.070934 

-191.063888 

-158.947068 

-158.946146 

-197.348647 

-197.351088 

Table  10:  Conformational  energies  in  atomic  units  for  three  representative  systems,  using  GVB, 
GVB-DFT  with  a  post-SCF  DFT  treatment,  and  our  preliminary  self-consistent  GVB-DFT 
software. 


Method 

formic  acid 

1-butene 

isoprene 

GVB 

4.31 

0.88 

1.34 

Post-SCF  DFT 

4.43 

0.57 

1.54 

SCF  GVB-DFT 

4.42 

0.58 

1.53 

Exp,  Value 

3.90 

0.53 

2.65 

Table  11:  Conformational  energy  differences  in  kcal/mol  for  three  representative  systems,  using 
GVB,  GVB-DFT  with  a  post-SCF  DFT  treatment,  and  our  preliminary  self-consistent  GVB- 
DFT  software. 


of  this  code  becomes  available,  we  will  be  able  to  calculate  atomization  and  bond  dissociation 
energies,  in  which  the  impact  of  a  fully  self-consistent  treatment  should  be  greater  (for 
instance,  the  RMS  error  for  the  atomization  energies  of  the  G2  database  molecules  is  halved 
for  the  standard  B3LYP  functional  when  it  is  applied  self-consistently  versus  as  a  post-SCF 
method). 

These  preliminary  results  are  encouraging  because  they  suggest  that,  for  at  least  some 
molecular  properties,  almost  all  of  the  accuracy  of  the  GVB-DFT  theory  can  be  obtained  by 
an  inexpensive  post-SCF  DFT  treatment.  Also,  when  a  self-consistent  calculation  is  required, 
as  for  the  calculation  of  analytic  gradients,  it  should  be  possible  to  devise  a  more  efficient  self- 
consistent  procedure  by  incorporating  the  DFT  exchange-correlation  potential  in  the  SCF 
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equations  only  after  the  orbitals  have  been  largely  converged  by  the  faster  standard  GVB 
method. 

In  Phase  II,  we  will  complete  the  development  of  the  SCF-GVB-DFT  software 
described  here  by  extending  it  to  open-shell  systems  and  mixed  “exact  exchange’VDFT 
methods.  This  will  be  the  basis  for  implementing  analytic  gradients  for  the  GVB-DFT  wave- 
function,  allowing  optimization  of  equilibrium  geometries  and  transition  states  at  this  level  of 
theory.  Existing  molecular  properties  methods  currently  available  for  standard  GVB  wave- 
functions  will  be  extended  to  GVB-DFT  wavefunctions  as  well. 
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